{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cbf709014ce0316e",
   "metadata": {},
   "source": [
    "## 基于单树GP的符号回归（Symbolic Regression）\n",
    "\n",
    "基于单树GP的符号回归是指使用遗传编程（GP）生成数学公式来逼近一组数据的关系，通过组合DEAP的Creator，Toolbox和Algorithms这三个模块即可实现。\n",
    "\n",
    "\n",
    "### Creator类\n",
    "Creator是一个工具类，其主要作用是创建新的类。在遗传编程中，通常需要自定义个体（Individual）和适应度（Fitness）类，因为不同的问题可能需要不同的适应度类型和个体结构。在DEAP中，我们可以使用creator来动态地创建这些类。\n",
    "\n",
    "在下面的例子中，我们创建了一个最基本的单目标单树GP，可以使用base.Fitness和gp.PrimitiveTree来定义。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "59cfefc0467c74ad",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-08T02:39:00.130308400Z",
     "start_time": "2023-11-08T02:39:00.012636500Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\zhenl\\anaconda3\\Lib\\site-packages\\deap\\creator.py:185: RuntimeWarning: A class named 'FitnessMin' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.\n",
      "  warnings.warn(\"A class named '{0}' has already been created and it \"\n",
      "C:\\Users\\zhenl\\anaconda3\\Lib\\site-packages\\deap\\creator.py:185: RuntimeWarning: A class named 'Individual' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.\n",
      "  warnings.warn(\"A class named '{0}' has already been created and it \"\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import operator\n",
    "\n",
    "from deap import base, creator, tools, gp\n",
    "\n",
    "\n",
    "# 符号回归\n",
    "def evalSymbReg(individual, pset):\n",
    "    # 编译GP树为函数\n",
    "    func = gp.compile(expr=individual, pset=pset)\n",
    "    # 计算均方误差（Mean Square Error，MSE）\n",
    "    mse = ((func(x) - x**2)**2 for x in range(-10, 10))\n",
    "    return (math.fsum(mse),)\n",
    "\n",
    "# 创建个体和适应度函数\n",
    "creator.create(\"FitnessMin\", base.Fitness, weights=(-1.0,))\n",
    "creator.create(\"Individual\", gp.PrimitiveTree, fitness=creator.FitnessMin)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "956e01e17271daa6",
   "metadata": {},
   "source": [
    "### Toolbox类\n",
    "Toolbox的作用类似于一个调度中心，它负责“注册”各种操作和函数。在遗传编程中，这些操作通常包括交叉（crossover）、变异（mutation）、选择（selection）和评估（evaluation）。通过register，我们可以将这些操作和相关的函数绑定在一起，以供后续算法使用。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "851794d4d36e3681",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-08T02:39:00.214209Z",
     "start_time": "2023-11-08T02:39:00.052073500Z"
    }
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "# 定义函数集合和终端集合\n",
    "pset = gp.PrimitiveSet(\"MAIN\", arity=1)\n",
    "pset.addPrimitive(operator.add, 2)\n",
    "pset.addPrimitive(operator.sub, 2)\n",
    "pset.addPrimitive(operator.mul, 2)\n",
    "pset.addPrimitive(operator.neg, 1)\n",
    "pset.addEphemeralConstant(\"rand101\", lambda: random.randint(-1, 1))\n",
    "pset.renameArguments(ARG0='x')\n",
    "\n",
    "# 定义遗传编程操作\n",
    "toolbox = base.Toolbox()\n",
    "toolbox.register(\"expr\", gp.genHalfAndHalf, pset=pset, min_=1, max_=2)\n",
    "toolbox.register(\"individual\", tools.initIterate, creator.Individual, toolbox.expr)\n",
    "toolbox.register(\"population\", tools.initRepeat, list, toolbox.individual)\n",
    "toolbox.register(\"compile\", gp.compile, pset=pset)\n",
    "toolbox.register(\"evaluate\", evalSymbReg, pset=pset)\n",
    "toolbox.register(\"select\", tools.selTournament, tournsize=3)\n",
    "toolbox.register(\"mate\", gp.cxOnePoint)\n",
    "toolbox.register(\"mutate\", gp.mutUniform, expr=toolbox.expr, pset=pset)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62f30d17704db709",
   "metadata": {},
   "source": [
    "### Algorithms类\n",
    "Algorithms模块提供了一些现成的遗传算法和遗传编程的实现。例如，eaSimple是一个简单的遗传算法，它可以处理基本的选择、交叉、变异和演化迭代。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "515b587d4f8876ea",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-08T02:39:00.216839200Z",
     "start_time": "2023-11-08T02:39:00.068850700Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   \t      \t                       fitness                        \t                      size                     \n",
      "   \t      \t------------------------------------------------------\t-----------------------------------------------\n",
      "gen\tnevals\tavg    \tgen\tmax        \tmin\tnevals\tstd   \tavg \tgen\tmax\tmin\tnevals\tstd    \n",
      "0  \t100   \t71177.5\t0  \t3.19748e+06\t0  \t100   \t314314\t4.04\t0  \t7  \t2  \t100   \t1.60574\n",
      "1  \t93    \t43553.5\t1  \t162664     \t0  \t93    \t29361.4\t4.65\t1  \t12 \t2  \t93    \t2.07545\n",
      "2  \t91    \t96589.4\t2  \t2.95681e+06\t0  \t91    \t398621 \t5.45\t2  \t12 \t2  \t91    \t3.01786\n",
      "3  \t94    \t2.53799e+06\t3  \t2.41098e+08\t0  \t94    \t2.39818e+07\t5.75\t3  \t17 \t2  \t94    \t3.00458\n",
      "4  \t91    \t199260     \t4  \t3.51947e+06\t0  \t91    \t704399     \t5.77\t4  \t14 \t2  \t91    \t2.85956\n",
      "5  \t92    \t2.48158e+06\t5  \t2.2959e+08 \t0  \t92    \t2.28354e+07\t5.57\t5  \t15 \t2  \t92    \t2.95044\n",
      "6  \t96    \t2.57693e+06\t6  \t2.35463e+08\t0  \t96    \t2.34184e+07\t5.93\t6  \t17 \t2  \t96    \t3.57563\n",
      "7  \t83    \t368231     \t7  \t3.51947e+06\t0  \t83    \t979332     \t5.67\t7  \t16 \t2  \t83    \t3.36171\n",
      "8  \t90    \t1.20119e+07\t8  \t9.41851e+08\t0  \t90    \t9.63415e+07\t5.54\t8  \t23 \t2  \t90    \t3.85855\n",
      "9  \t90    \t397038     \t9  \t3.19748e+06\t0  \t90    \t1.02046e+06\t5.37\t9  \t27 \t3  \t90    \t3.72734\n",
      "10 \t92    \t2.40345e+06\t10 \t2.2959e+08 \t0  \t92    \t2.28396e+07\t4.33\t10 \t14 \t3  \t92    \t2.14967\n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "from deap import algorithms\n",
    "\n",
    "# 定义统计指标\n",
    "stats_fit = tools.Statistics(lambda ind: ind.fitness.values)\n",
    "stats_size = tools.Statistics(len)\n",
    "mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size)\n",
    "mstats.register(\"avg\", numpy.mean)\n",
    "mstats.register(\"std\", numpy.std)\n",
    "mstats.register(\"min\", numpy.min)\n",
    "mstats.register(\"max\", numpy.max)\n",
    "\n",
    "# 使用默认算法\n",
    "population = toolbox.population(n=100)\n",
    "hof = tools.HallOfFame(1)\n",
    "pop, log  = algorithms.eaSimple(population=population,\n",
    "                           toolbox=toolbox, cxpb=0.9, mutpb=0.1, ngen=10, stats=mstats, halloffame=hof, verbose=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "237b39454ea988bc",
   "metadata": {},
   "source": [
    "由于DEAP重载了字符串运算符，因此可以直接输出结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "918142f4e60d65a0",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-08T02:39:00.217794500Z",
     "start_time": "2023-11-08T02:39:00.118939200Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mul(neg(x), neg(x))\n"
     ]
    }
   ],
   "source": [
    "print(str(hof[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54fe3d72a677307c",
   "metadata": {},
   "source": [
    "当然，我们也可以利用NetworkX库来对GP树进行可视化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "2fa44e7277d90c4c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-08T02:39:00.449935300Z",
     "start_time": "2023-11-08T02:39:00.134624200Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGFCAYAAABg2vAPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAz4ElEQVR4nO3de1hUdeLH8c8woqgBmhfUzdVtN7NdL8NFNK8UJZplqampGSnesbK8ZdZ2s6vZxUuihqKl5qXM1E0MQ00JhYHxsllrbrr+Mi8ZgqIQMvP7w9VdSldQhjMz5/16Hp8nhsP5fvAE34/ne+Yci8vlcgkAAJiWn9EBAACAsSgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMLlKpdnI6XTq8OHDCgwMlMVicXcmAABQDlwul06dOqUGDRrIz+/y//4vVRk4fPiwGjZsWG7hAABAxTl06JBuuOGGy36+VGUgMDDw4s6CgoLKJxkAAHCrvLw8NWzY8OI8fjmlKgMXlgaCgoIoAwAAeJkrLfFzASEAACZHGQAAwOQoAwAAmBxlAAAAk6MMAABgcpQBAABMjjIAAIDJUQYAADA5ygAAACZHGQAAwOQoAwAAmBxlAAAAk6MMAABgcqV6aiEA31FUVKxdu47Lbj+qrKxj+vHH0yosLFaVKlbVr3+dwsLqKjw8RC1a1JG/v9XouAAqAGUAMImDB3M1Z84uJSTsVE5OgSTJ399PRUXOi9v4+/tpzpzzH9esGaARI1pq+PAWatQo2JDMACqGxeVyua60UV5enoKDg5Wbm6ugoKCKyAWgnOTmFmrcuE1KTNwtPz+Liouv+CN/kdVqkdPpUlxcc02bFqWgoCpuTAqgvJV2/uaaAcCHbdhwQE2bzteCBXvkcqlMRUA6v73LJS1YsEdNm87Xhg0H3BMUgKEoA4CPmjkzSzExK3Xs2Jkyl4BfKy526ejRM4qJWalZs7LLKSEAT0EZAHzQrFnZeuSRLyRJTue1FYELLuxn9OiNFALAx1AGAB+zYcMBjR690a1jjB69kSUDwIdQBgAfkptbqNjYz+TnZ3HrOH5+Fj388GfKyyt06zgAKgZlAPAh48Zt0vHjZ8ptaeBynE6Xjh07o7FjN7l1HAAVgzIA+IgDB3KVmLj7mi8WLK3iYpcSE3fr4MHcChkPgPtQBgAfMXfuLvn5WRQb+xe5XOPkco1Tp04NL7ntvn1xcrnGKTW1b5nHefbZtnK5xkk6v1wwd+6ua8oNwHiUAcAHFBUVKyFhZ4mzAnl5hYqLa/abbTt1aqg//almuaz3Fxe7NHv2ThUVFV/zvgAYhzIA+IBdu45fvMXwBcuWfatevZooMLByidfj4popLe0H/etfp8pl7JycAu3e/VO57AuAMSgDgA+w24/+5rWlS/dKkvr1a3rxtaCgyurVq4nmz99TYttOnRpeclmhUaMguVzjFBv7lzKPD8B7UAYAH5CVdUz+/iV/nPPyftHKlf/Q4MHNL77Wr98tcjpdWrbsm3Ib29/fjzIAeDnKAOADfvzxdImnD14wf/4etW5dX3/+cy1J0uDBzbRixT90+nRRuY1dVOTUkSP55bY/ABWPMgD4gMLCS1/At3nzIX33XY4GD26uZs1qKzKyvubP313u4xcUnCv3fQKoOJWMDgDg2lWpYr3s5xYs2KNHHw1TQIBV3377s7Zu/aHcxw8I4FcJ4M04MwD4gPr1r/vNNQMXJCX9XbVrV9WIES21YMGeS25z4V/2vy4VtWtXveLY/v5+qlevehkTA/AklAHAB4SF1b3kNQOSdPjwaU2dmqE1a/Zr4cK/X3KbAwfO30WwRYs6JV7v3v1PVxy7qMip8PCQMiYG4Ek4twf4gCtNxpMmffk/P3/06Bl9/vkBTZoUqZycAh08mKfo6N+rZ8+bymV8AJ6NMwOAD2jRoo5q1gy4pn0MHPg3bdz4L732WketWHGPfve769Sv37orfl3NmgFq3rz2NY0NwFgWl8t1xaea5OXlKTg4WLm5uQoKCqqIXADK6KmnvtTrr++osAcVSZLVatHEiZF66aUOFTYmgNIr7fzNmQHARwwf3sLtjy7+NafTpWHDWlTomADKH2UA8BGNGgUrLq65/Crop9pqtSgurrkaNQqumAEBuA1lAPAhbdv+LJfrlKRLv7Og/DgVGGjRtGlRbh4HQEWgDAA+4OzZsxoxYoQGD+6vjh1/lPt/tP108uQcTZw4RmfPnnXzWADcjTIAeLlvv/1Wbdq00cKFCzVnzhylps7RzJnRbh1z5sxozZkzVgsWLFCbNm307bffunU8AO5FGQC82AcffKDw8HAVFhZq+/btGjZsmCwWi+LjQy8WAj8/S7mMdWE/s2ZFKz4+VMOGDdOOHTtUUFCg8PBwLV68uFzGAVDxKAOAFzpz5ozi4uI0cOBA9ezZU5mZmWrRouRV/fHxoUpOvl8hIdVktV5bIbBaLQoJqabk5Ps1alToxddbtGghu92uHj166MEHH9SQIUN05syZaxoLQMWjDABe5uuvv1arVq20dOlSzZ8/XwsXLtR11113yW07d26svXsHa9CgZrJYVOZSYLVaZLFIgwY10zffDFbnzo1/s811112nRYsWKTExUUuWLFFkZKS+/vrrq/nWABiEMgB4CZfLpQULFigiIkKSlJGRoUGDBsli+d8TfHBwFc2bF6Pvvx+qiRMjS9yp8NcPN/rvj2vWDNDEiZH6/vuhmjcvRkFBVS47hsVi0eDBg5WRkSGXy6VWrVopKSnpKr5LAEbgDoSAFzh9+rRGjRql999/X4MHD9aMGTNUrVq1q9pXUVGxdu/+SXb7UdntR3XkSL4KCs4pIKCS6tWrrvDwEIWHh6h589ry97/8o5EvJz8/X4888ogWLFighx56SLNmzbrsmQsA7lXa+ZsyAHi4Xbt2qW/fvjp06JASEhL04IMPGh2pVN5//32NHDlSDRs21PLly9W8eXOjIwGmw+2IAS/ncrk0d+5ctW7dWpUrV1ZmZqbXFAFJGjhwoDIzM+Xv76/IyEjNmzdPpfi3BwADUAYAD5SXl6f+/ftr+PDhio2NVXp6upo2bWp0rDJr2rSptm/froceekjDhg3TgAEDdOrUKaNjAfgVygDgYbKzsxUeHq5169Zp6dKlSkhIUNWqVY2OddWqVq2qOXPmaOnSpVqzZo3CwsKUnZ1tdCwA/4UyAHgIl8uld999V23atFFgYKCysrL0wAMPGB2r3DzwwAPKyspSYGCgbr31Vr377rssGwAegjIAeICTJ0+qT58+io+P17Bhw5SWlqY//elPRscqdzfddJPS0tI0ZMgQxcfHq0+fPsrNzTU6FmB6lAHAYBkZGQoLC9Pnn3+ulStXasaMGQoICLjyF3qpgIAAzZw5UytWrNCGDRsUGhqqzMxMo2MBpkYZAAzicrn09ttvq127dqpdu7aysrLUq1cvo2NVmPvvv1/Z2dmqVauW2rZtq3feeYdlA8AglAHAAD///LN69Oihxx9/XPHx8dq6datuvPFGo2NVuBtvvFHbtm1TfHy8xowZo549eyonJ8foWIDpUAaACpaenq7Q0FBt2bJFq1ev1ltvvaXKlSsbHcswlStX1ltvvaVPPvlEmzZtUmhoqNLT042OBZgKZQCoIE6nU1OnTlWHDh30u9/9TtnZ2erevbvRsTzGvffeK4fDofr166tDhw5644035HQ6jY4FmAJlAKgAP/30k7p3764JEyboiSee0ObNm9WoUSOjY3mcRo0aacuWLXr88cc1fvx4de/eXSdOnDA6FuDzKAOAm23dulU2m03p6elat26dXnvtNfn7+xsdy2P5+/vr9ddf19q1a5Weni6bzaatW7caHQvwaZQBwE2cTqdeeeUVRUVF6Q9/+IMcDofuuusuo2N5jW7dusnhcKhx48aKiorSK6+8wrIB4CaUAcANjh07pq5du2ry5Ml68sknlZqaqhtuuMHoWF7nhhtuUGpqqiZOnKjJkyera9euOnbsmNGxAJ9DGQDK2aZNm2Sz2ZSdna3k5GRNmTJFlSpVMjqW16pUqZJeeuklrV+/XtnZ2bLZbNq0aZPRsQCfQhkAyklxcbGef/55RUdHq2nTptq5c6fuvPNOo2P5jM6dO8vhcOjmm29WdHS0XnjhBRUXFxsdC/AJlAGgHBw5ckSdO3fW888/r2eeeUaff/656tevb3Qsn9OgQQOlpKTomWee0XPPPafOnTvryJEjRscCvB5lALhGKSkpatmypb7++mulpKToueeek9VqNTqWz7JarXruueeUkpKir7/+WjabTSkpKUbHArwaZQC4SufOndPTTz+tzp07q2XLlnI4HLr99tuNjmUat99+uxwOh5o3b67OnTvrmWee0blz54yOBXglygBwFX744QdFR0frlVde0ZQpU7R+/XqFhIQYHct0QkJClJycrBdffFEvv/yyoqOj9cMPPxgdC/A6lAGgjNavXy+bzab9+/dr06ZNeuqpp+Tnx4+SUfz8/DR58mSlpqbqu+++k81m0/r1642OBXgVfoMBpVRUVKQnn3xSXbt2VatWreRwONShQwejY+HfOnbsKIfDoYiICHXt2lWTJk1SUVGR0bEAr0AZAErhX//6l6KiovTGG29cvFVu7dq1jY6FX6lTp87FWz5PnTpVUVFROnTokNGxAI9HGQCuYM2aNQoNDdWhQ4f05Zdfavz48SwLeDA/Pz9NmDBBW7Zs0aFDh2Sz2bR27VqjYwEejd9owGX88ssvGjt2rLp376727dvL4XDo1ltvNToWSqlt27bKzs5Wu3btdM8992js2LH65ZdfjI4FeCTKAHAJ33//vTp06KAZM2bozTff1CeffKLrr7/e6Fgoo1q1amn16tV68803NX36dHXo0EEHDhwwOhbgcSgDwK98/PHHCg0N1bFjx7R161Y9/vjjslgsRsfCVbJYLHr88ce1bds2HTt2TKGhoVq1apXRsQCPQhkA/q2wsFCPPPKIevXqpejoaGVnZysyMtLoWCgnkZGRys7O1m233aaePXvq0UcfVWFhodGxAI9AGQAkfffdd2rbtq3mzp2rGTNmaOXKlapRo4bRsVDOatSooY8++kgzZszQnDlz1K5dO+3fv9/oWIDhKAMwveXLlyssLEy5ubn66quvNHr0aJYFfJjFYtHo0aOVlpamkydPKiwsTCtWrDA6FmAoygBM6+zZsxo5cqT69u2ru+66S1lZWQoLCzM6FipIeHi4srKy1KVLF/Xp00ejRo1SQUGB0bEAQ1AGYErffvut2rRpowULFmjOnDlaunSpgoKCjI6FChYUFKQPP/xQCQkJmj9/vtq0aaN//OMfRscCKhxlAKazePFihYeHq6CgQDt27NCwYcNYFjAxi8Wi4cOHa/v27Tp79qzCw8O1ZMkSo2MBFYoyANM4c+aMhgwZogcffFA9evSQ3W5XixYtjI4FD9GyZUtlZmbq3nvv1YABAzR06FCdOXPG6FhAhaAMwBS+/vprRUZGasmSJUpMTNSiRYt03XXXGR0LHiYwMFDvv/++EhMTtXjxYrVu3Vp79+41OhbgdpQB+LykpCS1atVKLpdLGRkZGjx4MMsCuCyLxaLBgwdrx44dKi4uVkREhBYuXGh0LMCtKAPwWadPn1ZsbKwGDRqkvn37aseOHfrLX/5idCx4iWbNmikjI0N9+vTRww8/rNjYWOXn5xsdC3ALygB80u7du9WqVSt99NFHWrRokebPn6/q1asbHQtepnr16lqwYIEWLlyolStXKiIiQrt37zY6FlDuKAPwKS6XS/PmzVNkZKT8/f2VmZmpgQMHGh0LXu6hhx6S3W5XpUqVFBkZqffee08ul8voWEC5oQzAZ5w6dUoDBgzQsGHD9NBDD2n79u1q2rSp0bHgI5o2baodO3Zo4MCBGjp0qB588EGdOnXK6FhAuaAMwCdkZ2crLCxMa9as0dKlSzVnzhxVrVrV6FjwMVWrVtXcuXO1ZMkSffrppwoPD5fD4TA6FnDNKAPwai6XS++++65uvfVWBQYGKisrSw888IDRseDj+vXrp6ysLFWvXl1t2rTR7NmzWTaAV6MMwGvl5uaqT58+io+P15AhQ5SWlqabbrrJ6FgwiZtuuklfffWV4uLiNGrUKPXt21e5ublGxwKuCmUAXikzM1NhYWHasGGDVqxYoZkzZyogIMDoWDCZgIAAzZo1SytWrFBycrLCwsKUmZlpdCygzCgD8Coul0vvvPOO2rZtq+uvv17Z2dm6//77jY4Fk7v//vuVnZ2t66+/Xm3bttX06dNZNoBXoQzAa+Tk5Khnz54aM2aM4uPjtW3bNt14441GxwIkSTfeeKO2bt2qUaNG6bHHHlOvXr2Uk5NjdCygVCgD8Arp6ekKDQ3Vpk2b9Mknn+itt95S5cqVjY4FlFClShW9/fbbWrVqlVJTUxUaGqrt27cbHQu4IsoAPJrT6dQbb7yhDh06qH79+nI4HLr33nuNjgX8T/fdd58cDofq1aun9u3ba9q0aSwbwKNRBuCxTpw4oe7du2v8+PF6/PHHtWXLFjVq1MjoWECpNGrUSF9++aXGjBmjcePGqXv37jpx4oTRsYBLogzAI23dulU2m03p6elau3atXn/9dfn7+xsdCygTf39/TZ06VWvXrtVXX30lm82mbdu2GR0L+A3KADyK0+nUK6+8oqioKDVu3FgOh0PdunUzOhZwTbp16yaHw6FGjRqpU6dOevXVV+V0Oo2OBVxEGYDHOHbsmO666y5NnjxZEydOVGpqqm644QajYwHl4oYbbtCmTZs0YcIETZo0Sd26ddPx48eNjgVIogzAQ2zevFk2m01ZWVlav369XnrpJVWqVMnoWEC5qlSpkl5++WWtX79edrtdNptNmzdvNjoWQBmAsYqLi/XCCy/o9ttv18033yyHw6HOnTsbHQtwq5iYGDkcDjVp0kS33367XnzxRRUXFxsdCyZGGYBhjhw5opiYGD333HN65plnlJKSogYNGhgdC6gQDRo0UEpKip5++mk9++yziomJ0ZEjR4yOBZOiDMAQGzdulM1m09///nelpKToueeek9VqNToWUKGsVquef/55ff7559qzZ49sNps2btxodCyYEGUAFercuXP661//qjvvvFPNmzeXw+HQ7bffbnQswFDR0dFyOBxq1qyZ7rzzTj377LMsG6BCUQZQYX744QdFR0frpZde0osvvqj169crJCTE6FiAR6hXr56Sk5P1wgsvaMqUKYqOjtbhw4eNjgWToAygQqxfv142m03fffedUlNTNXnyZJYFgF+xWq16+umn9cUXX2jfvn1q2bKlkpOTjY4FE6AMwK2Kioo0adIkde3aVREREXI4HOrYsaPRsQCP1qlTJzkcDkVERKhLly6aNGmSzp07Z3Qs+DDKANzm0KFDioqK0tSpU/Xaa69p3bp1qlOnjtGxAK9Qp04drVu3Tq+++qqmTp2qqKgoHTp0yOhY8FGUAbjF2rVrZbPZdOjQIW3ZskUTJkyQnx//uwFl4efnp4kTJ2rz5s06ePCgbDab1q1bZ3Qs+CB+O6Nc/fLLLxo7dqzuuecetWvXTtnZ2Wrbtq3RsQCv1q5dOzkcDrVt21Z33323xo0bp6KiIqNjwYdQBlBuDhw4oI4dO2r69Ol68803tXr1atWqVcvoWIBPqFWrlj799FNNmzZN77zzjjp06KADBw4YHQs+gjKAcvHJJ58oNDRUR48e1bZt2/T444/LYrEYHQvwKRaLRU888YS2bt2qI0eOKDQ0VJ988onRseADKAO4JoWFhXrsscfUo0cP3XbbbcrOzlZkZKTRsQCf1rp1a2VnZysqKko9evTQmDFjVFhYaHQseDHKAK7a/v371a5dOyUkJGjGjBn66KOPVKNGDaNjAaZQs2ZNffzxx5o+fbpmz56tdu3a6Z///KfRseClKAO4KitWrFBYWJhOnjyptLQ0jR49mmUBoIJZLBY98sgjSktLU05OjkJDQ7Vy5UqjY8ELUQZQJgUFBRo1apT69OmjLl26yG63Kzw83OhYgKmFh4crKytLMTEx6t27t0aNGqWCggKjY8GLUAZQav/4xz/Upk0bzZ8/X7Nnz9aHH36o4OBgo2MBkBQcHKxly5Zp9uzZmj9/vm699Vbt27fP6FjwEpQBlMqSJUsUHh6us2fPavv27RoxYgTLAoCHsVgsGjFihNLT05Wfn6+wsDAtXbrU6FjwApQB/E9nzpzR0KFDNWDAAN17773KzMxUy5YtjY4F4H+w2Wyy2+3q3r27+vfvr6FDh+rs2bNGx4IHowzgsvbu3avWrVtr8eLFSkxM1Pvvv6/AwECjYwEohcDAQH3wwQd677339MEHHygyMlJ79+41OhY8FGUAl7Rw4UJFRESouLhYO3bs0ODBg1kWALyMxWJRXFycMjIyVFxcrIiICC1atMjoWPBAlAGUkJ+fr4cfflgPP/yw+vTpo4yMDDVr1szoWACuQbNmzZSRkaHevXsrNjZWgwYNUn5+vtGx4EEoA7hoz549ioiI0IoVK7Rw4UItWLBA1atXNzoWgHJQvXp1JSUlKSkpScuXL1erVq20Z88eo2PBQ1AGIJfLpffee0+tWrVSpUqVZLfb9dBDDxkdC4AbxMbGKjMzU1arVZGRkUpMTJTL5TI6FgxGGTC5U6dO6cEHH9TQoUM1cOBA7dixQ02bNjU6FgA3uuWWW7R9+3YNGDBAQ4YM0cCBA3Xq1CmjY8FAlAETczgcioiI0KeffqolS5Zo7ty5qlq1qtGxAFSAatWqad68eVq8eLFWr16tiIgI7dy50+hYMAhlwIRcLpdmz56tNm3aqFq1arLb7erXr5/RsQAYoH///rLb7apatapat26thIQElg1MiDJgMrm5uerbt69GjRqluLg4ffXVV2rSpInRsQAYqEmTJkpPT1dcXJxGjhypBx54QLm5uUbHQgWiDJhIZmamwsLClJycrBUrVmjWrFkKCAgwOhYADxAQEKBZs2Zp+fLlWr9+vcLDw2W3242OhQpCGTABl8ul6dOnq23btrr++uuVnZ2t+++/3+hYADxQ7969lZWVpRo1aqht27aaMWMGywYmQBnwcTk5OerVq5cee+wxjRo1Slu3btWNN95odCwAHuyPf/yjtm3bphEjRujRRx9Vr169lJOTY3QsuBFlwIdt375doaGhSk1N1apVq/T222+rSpUqRscC4AWqVKmid955R6tWrVJqaqrCwsK0Y8cOo2PBTSgDPsjlcmnatGlq37696tWrJ4fDofvuu8/oWAC80H333afs7GyFhISoXbt2evPNN1k28EGUAR9z4sQJde/eXePGjdOYMWP05ZdfqlGjRkbHAuDFGjdurC1btuixxx7T2LFjde+99+rnn382OhbKEWXAh2zbtk2hoaFKS0vT2rVrNXXqVPn7+xsdC4APqFy5st544w2tWbNG27Ztk81mU1pamtGxUE4oAz7A6XTq1VdfVadOnfT73/9eDodD3bp1MzoWAB909913y+Fw6Pe//706duyo1157TU6n0+hYuEaUAS93/PhxdevWTZMmTdKECRO0adMmNWzY0OhYAHxYw4YNlZqaqvHjx+vJJ5/U3XffrePHjxsdC9eAMuDFNm/eLJvNJrvdrvXr1+vll19WpUqVjI4FwAT8/f31yiuv6LPPPlNGRoZsNpu2bNlidCxcJcqAFyouLtaLL76o22+/XU2aNJHD4VBMTIzRsQCYUJcuXbRz507ddNNNuu222zRlyhQVFxcbHQtlRBnwMkeOHFFMTIyeffZZPf3000pJSVGDBg2MjgXAxBo0aKCUlBRNnjxZf/3rX9WlSxcdPXrU6FgoA8qAF9m4caNsNpv27Nmjzz//XM8//7ysVqvRsQBAlSpV0gsvvKANGzZo9+7datmypb744gujY6GUKANeoLi4WM8++6zuvPNONWvWTA6HQ9HR0UbHAoDfuOOOO+RwONSsWTPdcccdevbZZ1k28AKUAQ93+PBhRUdHa8qUKXrhhReUnJysevXqGR0LAC6rXr16Sk5O1vPPP68pU6bojjvu0OHDh42Ohf+BMuDBkpOTZbPZtG/fPn3xxRd6+umnWRYA4BWsVqueeeYZffHFF/r2229ls9m0YcMGo2PhMigDHujcuXOaNGmSunTpovDwcDkcDnXq1MnoWABQZp06dZLD4VBYWJhiYmL01FNP6dy5c0bHwq9QBjzMoUOHFBUVpalTp+rVV1/VunXrVKdOHaNjAcBVq1u3rv72t7/plVde0euvv67bbrtN//d//2d0LPwXyoAHWbdunWw2mw4ePKjNmzdr4sSJ8vPjEAHwfn5+fnryySe1adMmHThwQDabTX/729+MjoV/Y6bxAEVFRRo/frzuvvtutW3bVg6HQ+3atTM6FgCUu/bt28vhcOjWW29Vt27dNGHCBBUVFRkdy/QoAwY7cOCAOnTooLffflvTpk3Tp59+qlq1ahkdCwDcplatWvr000/1xhtv6K233lLHjh118OBBo2OZGmXAQJ988olCQ0N15MgRbd26VU888YQsFovRsQDA7SwWi8aOHasvv/xSP/74o2w2m1avXm10LNOiDBigsLBQY8aMUY8ePRQVFaXs7Gy1bt3a6FgAUOHatGmj7OxsRUVF6b777tOYMWP0yy+/GB3LdCgDFeyf//yn2rVrp9mzZ2v69On6+OOPVbNmTaNjAYBhatasqY8//ljvvPOO3n33XbVr107//Oc/jY5lKpSBCrRy5UqFhoYqJydHaWlpeuSRR1gWAACdXzZ49NFHlZaWpp9//lmhoaFauXKl0bFMgzJQAQoKChQfH6/evXsrJiZGWVlZCg8PNzoWAHiciIgIZWVlqXPnzurdu7fi4+NVUFBgdCyfRxlws3379unWW29VYmKiZs+erWXLlik4ONjoWADgsYKDg7V8+XK9++67SkxMVNu2bbVv3z6jY/k0yoAbLV26VGFhYcrPz1d6erpGjBjBsgAAlILFYtHIkSOVnp6u06dPKywsTB9++KHRsXwWZcANzp49q2HDhql///7q3r277Ha7bDab0bEAwOvYbDbZ7Xbdc8896tevn4YPH66zZ88aHcvnUAbK2d69exUZGan3339f7733nj744AMFBgYaHQsAvFZgYKAWL16sefPmadGiRWrdurW++eYbo2P5FMpAOVq0aJEiIiJUXFysjIwMxcXFsSwAAOXAYrFoyJAh2rFjh4qKihQeHq5FixYZHctnUAbKQX5+vgYNGqTY2Fj17t1bGRkZatasmdGxAMDnNG/eXJmZmerdu7diY2M1aNAg5efnGx3L61EGrtGePXvUqlUrLV++XElJSUpKSlL16tWNjgUAPqt69eoXf98uX75ckZGR+vvf/250LK9GGbhKLpdLiYmJioyMlNVqVWZmpmJjY42OBQCmERsbq4yMDFksFrVq1Urz58+Xy+UyOpZXogxchVOnTmngwIEaMmSIBgwYoO3bt+uWW24xOhYAmM6f//xn7dixQwMGDFBcXJwGDhyo06dPGx3L61AGymjnzp2KiIjQ6tWrL17dWq1aNaNjAYBpVatWTfPmzdPixYu1evVqhYeHa+fOnUbH8iqUgVJyuVxKSEhQ69atVbVqVdntdvXv39/oWACAf+vfv7/sdruqVq2q1q1ba86cOSwblBJloBTy8vL0wAMPaOTIkRo8eLDS09PVpEkTo2MBAH6lSZMmSk9P16BBgzRixAj169dPeXl5RsfyeJWMDnCtioqKtWvXcdntR5WVdUw//nhahYXFqlLFqvr1r1NYWF2Fh4eoRYs68ve3lnn/drtdffv21fHjx7Vs2TL16dPHDd8FAKC8BAQEaPbs2brttts0ZMgQhYWFafny5QoLC7uq/bl7nvEEFlcpzqHk5eUpODhYubm5CgoKqohcV3TwYK7mzNmlhISdysk5/0Qrf38/FRU5L27z3x/XrBmgESNaavjwFmrU6MoPCnK5XJo5c6bGjRun5s2ba9myZfrjH//onm8GAOAW+/fvV58+fbRnzx5NmzZN8fHxpb4ZnLvnmYpQ2vnb68pAbm6hxo3bpMTE3fLzs6i4uPTrQVarRU6nS3FxzTVtWpSCgqpccrucnBzFxcVp1apVevTRR/X666+rSpVLbwsA8GyFhYUaP368ZsyYoZ49eyoxMVE1atS47PYVMc9UFJ8sAxs2HFBs7Gc6fvxMmQ7Or1mtFtWtW01JSV3VuXPjEp/bsWOH+vbtq5MnT2rBggW67777ri00AMAjrFq1SoMHD1aNGjW0bNkyRUZG/mabiphnKlJp52+vuYBw5swsxcSs1LFj13aAJKm42KWjR88oJmalZs3KlnR+WeDNN99Uu3btFBISouzsbIoAAPiQHj16KDs7W3Xr1lX79u311ltvlXi3gbvnGU/mFRcQzpqVrUce+UKS5HSWz9tELuxn9OiNys/P19atr2vNmjUaO3asXn75ZVWuXLlcxgEAeI7GjRvryy+/1FNPPaUnnnhCqampSkpK0tKlB906z0hSfHxouezXHTx+mWDDhgOKiVnp9nECAz/UkiXP6e6773b7WAAA461Zs0YPP/ywrNZbdPz4fW4fLzn5/gpfMvCJZYLc3ELFxn4mPz93PwbYqapVH1bHjne6eRwAgKe45557tHnzdp082UWS84rbXws/P4sefvgz5eUVunWcq+XRZWDcuE06fvxMuZ2yuTw/nThRqLFjN7l5HACAJ3nnnf1yOqvK3dOh0+nSsWNnPHae8dhlgp9+OqPataupoOCcbr55vv71r5J3kEpN7avatauqefOkchvTYpG+/36ox7w/FADgPmaYZ7x+mWDTpkOSpICASpoypV2FjOnnZ9HcubsqZCwAgLGYZ/7DI8tAUVGxUlPPH6TPPvte/fvfohYt6rh93OJil2bP3qmiomK3jwUAMA7zTEkeWQZ27Tqu/PwiSdLrr+/QiRMFeu21jlf8upEjbcrOfkhnzjymn38erRUruusPf/jtqZhJk1rrwIFhOnt2jDIyHtQddzRSampfpab2VU5OgXbv/qncvycAgOdgninJI8uA3X704n+fOvWLpkz5Sl26/EG33dbwsl8zZ86devvt25SSclD33bdao0al6C9/qaW0tP6qW7faxe1eeqm9Xn65g9av/1733vuJEhJ26r33YtSkSc1Ljg8A8D3MMyV5ZBnIyjomq/U/bydMSNip/ftP6rXXOl1y+9at62vYsJZ68sktGj9+szZsOKAPP/xGd965QkFBlfXEExGSpBo1quiJJyL04YffaMSIz7VhwwElJu5W375r1KDBdZLOP3TC0w4SAKB8Mc+U5JFl4McfT5e4FWRRkVNPP71VrVrVU58+N/9m+7vvvlFOp0sffLBXVqvl4p8jR/K1c+dxRUWdb3pt2jRQQEAlLV/+bYmv3779R33/fe7FsY4cyXfjdwcAMBrzTEkeeTviwsLfXljx4YffaNy4CL30Unt9/PG+Ep8LCakuPz+Ljh0bdcn97d9/UpJUq1aAJOno0d8ehP9+raDg3NVGBwB4AeaZkjyyDFSpYr3k6xMnblFKSh8NG9aixOs//XRWTqdLHTosveQBvvDaiRPnn0cdElL9N9vUq1ddBw6cf49pQIBH/rUAAMoJ80xJHrlMUL/+dSXWci7YuPFf2rDhgP7611t13XX+F19fu3a//Pws+t3vAmW3H/3Nnz17zl+1uX37jyooOKe+fUueAmrdur4aNz5/Nai/v5/q1fvtQQQA+A7mmZI8q5r8W1hYXaWlHb7k5yZO3CK7faBCQqpf/MtPSzusOXN2asGCLoqICNGWLf+n/Pwi1a9fXe3b36Ddu48rIWGncnIK9OabmXrqqTbKySnQqlXf6YYbrtOzz7bV4cOn5XS6VFTkVHh4SEV+uwCACsY8U5JHloHw8JDLHiSH45iWLt2rAQP+XOL1ESM+V3r6jxo+vIVGjbLJz8+iw4dPa9u2w9qx48jF7SZP3qr8/CKNGNFSgwY10zff/KyRI1P00kvtdfJkwcXxAQC+i3mmJI98NkFRUbFCQmYrJ6fA7WNJUuPGwfrmm0F6/vmvlJCwU0ePjpS//6XXkwAA3s8s80xp52+PPDPg72/ViBEt9frrO0q89aM8tGhRR/36NVVa2mHl5f2im2+uqQkTIpWX94uSkvZo5MiWFAEA8HHMMyV55JkBSTp4MFd/+MM8XTld2fzxjzWUkHCnWrasoxo1qig3t1CbNv2fJk/+Uvv25fDUQgAwCTPMM159ZkCSGjUKVlxccy1YsKdcW9v+/Sd1550rfvO61WpRXFxzigAAmATzzH945FsLL5g2LUp161aTn99v3/5Rnvz8LKpbt5qmTYty6zgAAM/CPHOeR5eBoKAqSkrqKqeznM/h/IrT6VJSUlcFBVVx6zgAAM/CPHOeR5cBSercubFmzox26xizZkWrc+fGbh0DAOCZmGe8oAxIUnx86MUDVV6nci7sZ9asaI0aFVou+wQAeCezzzNeUQak8wcqOfl+hYRUu+QtJMvCarUoJKSakpPv9/gDBACoGGaeZ7ymDEjnT+Xs3TtYgwY1k8WiMh8sq9Uii0X/viPUYI8+ZQMAqHhmnWc89j4DV3LwYK7mzt2l2bN3XryDlL+/n4qKnBe3+e+Pa9YM0MiRLTVsWAuPfFsHAMCz+MI8U9r522vLwAVFRcXavfuni0+O+uabH7R58zZ16tROTZv+TuHhIQoPD1Hz5rU97o5PAADP583zjGnKwK9lZWUpPDxcdrtdYWFhRscBAPgYb5pnSjt/e9U1AwAAoPxRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJkcZAADA5CgDAACYHGUAAACTowwAAGBylYwOcK2Kioq1a9dx2e1HlZV1TN9884OkIXriiV1q2vSEwsLqKjw8RC1a1JG/v9XouAAAL2OGecbicrlcV9ooLy9PwcHBys3NVVBQUEXkuqKDB3M1Z84uJSTsVE5OgSTJ399PRUXOi9v898c1awZoxIiWGj68hRo1CjYkMwDAe/jCPFPa+dvrykBubqHGjdukxMTd8vOzqLj4ivEvslotcjpdiotrrmnTohQUVMWNSQEA3siX5hmfLAMbNhxQbOxnOn78TJkOzq9ZrRbVrVtNSUld1blz4/ILCADwar42z5R2/vaaCwhnzsxSTMxKHTt2bQdIkoqLXTp69IxiYlZq1qzsckoIAPBmZp5nvKIMzJqVrUce+UKS5HRe2wG64MJ+Ro/e6BUHCgDgPmafZzy+DGzYcECjR2906xijR2/Uhg0H3DoGAMAzMc94eBnIzS1UbOxn8vOzuHUcPz+LHn74M+XlFbp1HACAZ2GeOc+jy8C4cZt0/PiZcjtlczlOp0vHjp3R2LGb3DoOAMCzMM+c57Fl4MCBXCUm7r7mizhKq7jYpcTE3Tp4MLdCxgMAGIt55j88tgzMnbvL7adtfs3Pz6K5c3dV6JgAAGMwz/yHR5aBoqJiJSTsrLC2dkFxsUuzZ+9UUVFxhY4LAKhYzDMleWQZ2LXr+MVbP5ZGlSpWZWUN1L59cQoKqnzx9ZCQavrxx5FKTe1b6vaXk1Og3bt/KnNmAID3YJ4pySPLgN1+tEzbFxYWq0+fNapbt5rmz+8iSbJYpMWLu8likfr1W1umi0PKOj4AwLswz5TkkU8tzMo69puHQVzJd9+d1JAhG7R8+T169NEwXX99gKKiGqpLl4905Eh+qffj7+8nu/2ohg69muQAAG/APFOSR5aBH388XaYDdMGKFd/q3XcbaurUTrJaLXr55e1KSTlYpn0UFTnLdFABAN6HeaYkj1wmKCy8+gsr5s/frcqVrTp3zqnp07Ouah8FBeeuenwAgOdjninJI8tAlSrWq/q6atX89f77d+nbb3/W2bPn9N57MVe1n4AAjzxhAgAoJ8wzJXlkGahf/zr5+5c9WkLCHfr974PUs+dqxcUl6957/6QxY8LLtA9/fz/Vq1e9zGMDALwH80xJHlkGwsLqlnktJy6uuQYO/Ivi41P09dcn9PHH+zRjRpZee62jWrWqV+r9FBU5FR4eUtbIAAAvwjxTkkeWgbL+JTVrVlvTp9+upKQ9Wrjw7xdfHzdus3btOq5ly+5RcHAVt40PAPAuzDMlWVwu1xXfGJmXl6fg4GDl5uYqKCjI7aGKiooVEjK7TDeEKC81awbo6NGR8ve/uvUkAIDnM8s8U9r52yPPDPj7WzViREtZrRV7z2ir1aKRI1tSBADAxzHPlOSRZUCShg9v4fZHSv6a0+nSsGEtKnRMAIAxmGf+w2PLQKNGwYqLa15hrc1qtSgurrkaNQqukPEAAMZinvkPjy0DkjRtWpTq1q3m9kdM+vlZVLduNU2bFuXWcQAAnoV55jyPLgNBQVWUlNTV7adxnE6XkpK6Kiio9FeCAgC8H/PMeR5dBiSpc+fGmjkz2q1jzJoVrc6dG7t1DACAZ2Ke8YIyIEnx8aEXD1R5ncq5sJ9Zs6I1alRouewTAOCdzD7PeEUZkM4fqOTk+xUSUu2aL/awWi0KCamm5OT7Pf4AAQAqhpnnGa8pA9L5Uzl79w7WoEHNZLGozAfLarXIYpEGDWqmb74Z7NGnbAAAFc+s84xH3oGwNA4ezNXcubs0e/bOi3eQ8vf3K3Gv6f/+uGbNAI0c2VLDhrXwyLd1AAA8iy/MM6Wdv722DFxQVFSs3bt/kt1+VHb7UR05kq+CgnMKCKikevWqKzw8ROHhIWrevLbH3fEJAOD5vHmeMU0ZAAAAl+bVzyYAAAAVhzIAAIDJUQYAADA5ygAAACZHGQAAwOQoAwAAmBxlAAAAk6MMAABgcpQBAABMjjIAAIDJUQYAADA5ygAAACZHGQAAwOQqlWajCw82zMvLc2sYAABQfi7M21d6QHGpysCpU6ckSQ0bNrzGWAAAoKKdOnVKwcHBl/28xXWluiDJ6XTq8OHDCgwMlMViKdeAAADAPVwul06dOqUGDRrIz+/yVwaUqgwAAADfxQWEAACYHGUAAACTowwAAGBylAEAAEyOMgAAgMlRBgAAMDnKAAAAJvf/xyVVZtpydCAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import networkx as nx\n",
    "from deap.gp import graph\n",
    "from networkx.drawing.nx_agraph import graphviz_layout\n",
    "\n",
    "function_name = {\n",
    "    'add':'Add',\n",
    "    'sub':'Sub',\n",
    "    'mul':'Mul',\n",
    "    'neg':'Neg'\n",
    "}\n",
    "\n",
    "def is_number(string):\n",
    "    try:\n",
    "        float(string)\n",
    "        return True\n",
    "    except ValueError:\n",
    "        return False\n",
    "\n",
    "\n",
    "def plot_a_tree(tree=hof[0]):\n",
    "    red_nodes = []\n",
    "    purple_nodes = []\n",
    "    blue_nodes = []\n",
    "    for gid, g in enumerate(tree):\n",
    "        if (\n",
    "                hasattr(g, \"value\")\n",
    "                and isinstance(g.value, str)\n",
    "                and g.value.startswith(\"ARG\")\n",
    "        ):\n",
    "            g.value = g.value.replace(\"ARG\", \"X\")\n",
    "\n",
    "        if g.name in function_name:\n",
    "            g.name = function_name[g.name]\n",
    "\n",
    "        if hasattr(g, \"value\") and (\n",
    "                is_number(g.value)\n",
    "                or (g.value.startswith(\"X\") and int(g.value[1:]) < X.shape[1])\n",
    "        ):\n",
    "            # 基础节点\n",
    "            red_nodes.append(gid)\n",
    "        elif hasattr(g, \"value\") and g.value.startswith(\"X\"):\n",
    "            g.value = \"$\\phi$\" + str(int(g.value.replace(\"X\", \"\")) - X.shape[1] + 1)\n",
    "            purple_nodes.append(gid)\n",
    "        elif hasattr(g, \"value\") and g.value.startswith(\"$\\phi$\"):\n",
    "            purple_nodes.append(gid)\n",
    "        else:\n",
    "            # 深蓝色节点\n",
    "            blue_nodes.append(gid)\n",
    "    nodes, edges, labels = graph(tree)\n",
    "    g = nx.Graph()\n",
    "    g.add_nodes_from(nodes)\n",
    "    g.add_edges_from(edges)\n",
    "    pos = graphviz_layout(g, prog=\"dot\")\n",
    "    red_nodes_idx = [nodes.index(n) for n in nodes if n in red_nodes]\n",
    "    purple_nodes_idx = [nodes.index(n) for n in nodes if n in purple_nodes]\n",
    "    blue_nodes_idx = [nodes.index(n) for n in nodes if n in blue_nodes]\n",
    "    nx.draw_networkx_nodes(\n",
    "        g, pos, nodelist=red_nodes_idx, node_color=\"darkred\", node_size=500\n",
    "    )\n",
    "    nx.draw_networkx_nodes(\n",
    "        g, pos, nodelist=purple_nodes_idx, node_color=\"indigo\", node_size=500\n",
    "    )\n",
    "    nx.draw_networkx_nodes(\n",
    "        g, pos, nodelist=blue_nodes_idx, node_color=\"darkblue\", node_size=500\n",
    "    )\n",
    "    nx.draw_networkx_edges(g, pos)\n",
    "    nx.draw_networkx_labels(g, pos, labels, font_color=\"white\")\n",
    "\n",
    "\n",
    "plot_a_tree(hof[0])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
